RNA-Seq Analysis Unraveling Novel Genes and Pathways Influencing Corneal Wound Healing

Purpose Transdifferentiation of corneal fibroblasts to myofibroblasts in the stroma is a central mechanistic event in corneal wound healing. This study sought to characterize genes and pathways influencing transdifferentiation of human corneal fibroblasts (hCSFs) to human corneal myofibroblasts (hCMFs) using RNA sequencing (RNA-seq) to develop comprehensive mechanistic information and identify newer targets for corneal fibrosis management. Methods Primary hCSFs were derived from donor human corneas. hCMFs were generated by treating primary hCSFs with transforming growth factor β1 (TGFβ1; 5 ng/mL) for 72 hours under serum-free conditions. RNA was extracted using the RNeasy Plus Mini Kit and subjected to RNA-seq analysis after quality control testing. Differential gene expression, pathway enrichment, and protein–protein network analyses were performed using DESeq2, GSEA/PANTHER/Reactome, and Cytoscape/cytoHubba, respectively. Results RNA-seq analysis of hCMFs and hCSFs identified 3843 differentially expressed genes and transcripts (adjusted P < 0.05). The log(fold change) ≥ ±1.5 filter showed 816 upregulated and 739 downregulated genes between two cell types. Pathway enrichment analysis showed the highest normalized enrichment score for epithelial-to-mesenchymal transition (5.569), followed by mTORC1 signaling (2.949), angiogenesis (2.176), and TGFβ signaling (2.008). Protein–protein interaction network analysis identified the top 20 nodes influencing corneal myofibroblast development. The expression of a novel MXRA5 in corneal stroma and its association with corneal fibrosis was verified by real-time quantitative reverse transcription PCR and immunofluorescence. RNA-seq and gene count files were submitted to the NCBI Gene Expression Omnibus (GSE260476). Conclusions This study identified several novel genes involved in myofibroblast development, offering potential targets for developing newer therapeutic strategies for corneal fibrosis.

T he cornea, an outermost transparent, avascular, and immune-privileged tissue of the eye, contributes twothirds of refraction and protection to eyes. 1 An uncontrolled corneal healing after trauma or injury often leads to corneal scarring or fibrosis and vision loss in millions of people worldwide. 2][5] The stroma constitutes 90% of the cornea and contains keratocytes, collagens, extracellular matrix (ECM), and proteoglycans. 6Keratocytes, the quiescent transparent cells, play an important role in corneal homeostasis and wound healing. 7Keratocytes transdifferentiate into myofibroblasts under the influence of cytokines and facilitate tissue repair primarily after trauma or injury.Transforming growth factor β1 (TGFβ1), among many cytokines, has been shown to play a major role in the generation of myofibroblasts, which are characterized by de novo alpha-smooth muscle actin (αSMA) and collagen (Col) I and Col III proteins.Corneal stromal myofibroblasts express high amounts of intermediate filaments, fibronectin, vimentin, and desmin at the early stages of transdifferentiation. 8 Inadequate reductions and the persistence of corneal myofibroblasts in stroma after wound closure leads to corneal scarring and fibrosis and abnormal stromal ECMs, [9][10][11] and reduced crystallin alters the refractive index and light-scattering property of the cornea. 12igh-throughput next-generation RNA sequencing (RNAseq) and bioinformatics analyses allow identification of differentially expressed genes (DEGs), hub genes, and signaling pathways associated with corneal scar develop-ment. 13This technique is particularly useful for unraveling molecular complexities and discovering prospective targets for corneal scar treatment.This study examined transcriptional differences between human corneal stromal fibroblasts (hCSFs) and TGFβ1-induced fibrotic human corneal myofibroblasts (hCMFs) using RNA-seq and bioinformatics analyses.

METHODS Primary Culture of Human Corneal Stromal Cells
Healthy donor human corneas (three males and one female, 49-76 years of age) obtained from Saving Sight (Kansas City, MO, USA) and normal (non-fibrotic) and fibrotic cadaver corneas from human subjects were used.The use of such corneas does not constitute human subjects research according to the Department of Health and Human Services regulatory definitions.Primary hCSFs were generated by rinsing corneas with serum-free Minimum Essential Medium (MEM), removing the corneal endothelium and epithelium with a surgical scalpel blade, cutting them into small fragments, and incubating at 37°C in a 5% CO 2 incubator for 3 to 5 weeks in MEM supplemented with 10% heat-inactivated fetal bovine serum (Thermo Fisher Scientific, Waltham, MA, USA).

Differentiation of Fibroblasts to Myofibroblasts
The corneal fibrosis, in vitro, was initiated by inducing transdifferentiation of hCSFs into hCMFs by growing cultures in serum-free conditions for 8 hours followed by TGFβ1 (5 ng/mL) (PeproTech, Cranbury, NJ, USA) exposure under serum-free MEM for 72 hours.Cultures were fed with fresh media and TGFβ1 every 24 hours.

Immunofluorescence
Double immunofluorescence was used to verify levels of fibroblast-specific protein 1 (FSP1; a fibroblast marker) and αSMA (a myofibroblast marker) in hCSFs and hCMFs using a FIGURE 1. Schematic representation of the methodology adopted in this study for RNA-seq data analysis.The raw sequencing fragment numbers (read pairs) are shown as raw pairs, and fragment numbers after 5 ,3 -adapter trimmed and filtered <20-bp reads are shown as trimmed.The percentage of reads number aligned and not aligned to reference genome in trimmed pairs are shown as mapped and unmapped, respectively.

RNA Extraction and RNA-Seq
Primary hCSF and hCMF cultures were used in quadruplicate.To obtain RNA, cultures were first washed with 1× PBS, scraped with a "rubber policeman" in lysis buffer and subjected to the RNeasy Mini Kit (QIAGEN, Hilden, Germany) following the manufacturer's instructions.Isolated total RNA was reverse transcribed to cDNA using a commercial kit (Promega Corporation, Madison WI, USA) following the manufacturer's instructions.The reaction mixture was heated at 42°C for 30 minutes followed by enzyme inactivation at 90°C for 2 minutes. 14The RNA samples were sequenced by Arraystar (Rockville, MD, USA), and cDNAs were used for real-time quantitative reverse transcription PCR (qRT-PCR).

Differential Expression Analysis
The quadruplet samples of hCSF and hCMF were subjected to differential gene expression (DGE) analysis using a published method 15 and employing Python and R packages (R Foundation for Statistical Computing, Vienna, Austria) and bioinformatics tools.Figure 1 shows the approach utilized in this study.FastQC 0.12.1 16 and Cutadapt 4.6 17 were used for quality checks and to filter out adapter sequences, primers, poly-A tails, and other types of unwanted sequences from the high-throughput sequencing reads.We evaluated the quality of all eight RNA sequences from the four hCSF samples and four hCMF samples, and unwanted sequences were removed to ensure high-quality reads during the quality check.HISAT2 18 was used for readalignment of trimmed RNA sequences of both the hCSF and hCMF cell types with the human reference genome GRCh37.HISAT2 utilizes a novel genome indexing scheme by implementing a graph-based approach to capture a wide representation of genetic variant reads.HISAT2 implements the Graph FM (GFM) index, unlike other graph aligner methods (e.g., vg19, bpa aligner20). 19StringTie 2.2.1 20 was used for assembling mapped reads into transcripts via a network flow algorithm and optional de novo assembly.StringTie creates comprehensive and accurate gene reconstructions and precise estimations of expression levels compared to other prominent transcript assembly tools, such as Cufflinks, 21 IsoLasso, 22 Scripture, 23 and Traph. 24Furthermore, for each sample, StringTie merge was run with all assembled transcript files (in GTF/GFF format), reference annotation files (gencode.v32.annotation.gtf),and a global and unified set of transcripts (isoforms) from various RNAseq samples.The gene count files obtained from StringTie were used as input for DESeq2 1.34.0 to identify differentially expressed features. 25In addition, normalized gene counts were generated for gene set enrichment analysis (GSEA; version 4.3.2). 26 To add annotation information from the GTF file to the DEGs, Annotate DESeq2/DEXSeq was used.

Pathway Enrichment Analysis
To distinguish biological processes, cellular components, molecular functions, protein classes, and pathways of DEGs, protein annotation evolutionary relationship (PANTHER 18.0) classification (https://pantherdb.org/)and the Reactome database (https://reactome.org/) were used. 27A significance level of P < 0.05 was established as the threshold for determining substantial enrichment.To identify key pathways that are enriched in hCMFs compared to hCSFs, GSEA is a useful computational method that determines FIGURE 3. Volcano plot illustrating differential gene expression in the hCMF and hCSF cell types.This plot highlights significantly altered genes, which were upregulated and downregulated (in red) at a significance level of P ≤ 0.05 and a cutoff value of log 2 (FC) ± 1.5.The non-significant genes are depicted in gray and green, and genes outside the log 2 (FC) cutoff values are shown in blue.The plot distinctly shows the differential expression patterns between the hCMF and hCSF cell types, emphasizing the genes that are notable outliers in their expression levels.predefined statistically significant and concordant differences between two biological states (hCMFs and hCSFs).The enriched pathways were organized based on their normalized enrichment scores, and only those with a significance level of P < 0.01 were selected for subsequent analysis.

Protein-Protein Interaction Networks
To analyze known proteins and protein-protein interaction (PPI) networks, the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database and webbased search engine were used.The direct and indirect connections between proteins and their functional correlations were analyzed (https://string-db.org/)using a reported method. 28The STRING database was used for creating the corresponding PPI networks from the most enriched pathways obtained from the GSEA.Hub genes were identified based on network topology using the PPI network.Afterward, using Cytoscape 3.8.2 and the cytoHubba plugin, we identified the probable targets and influential nodes inside intricate networks. 29The analysis focused on examining the crucial hub genes that may be involved in the process of corneal wound healing.

Validation of Selected Genes Using qRT-PCR
For qRT-PCR, we used the StepOnePlus Real-Time PCR System (Applied Biosystems, Waltham, MA, USA) and gene-specific primers (Table 1) to validate expression of selected genes.A 20-μL reaction mixture contained 1 μL cDNA, 1 μL forward primer, 1 μL reverse primer, and 10 μL iQ SYBR1 Green Supermix (Bio-Rad Laboratories, Hercules, CA, USA) per reaction well.The mixture was exposed to the following PCR parameters: 95°C for 5 minutes, followed by 40 cycles of 95°C for 15 seconds, then 60°C for 1 minute, and then a final cycle of 72°C for 10 minutes.The fluorescence threshold value (C t ) was calculated to detect signal differences in association with an exponential increase of PCR products in the log-linear phase.The GAPDH gene was used for the normalization of data.Relative expression/fold change (FC) over the corresponding values for the control was calculated by the 2 -Ct method.Two or three independent experiments were executed; for each sample, qRT-PCR was performed in triplicate, and the average fold changes in mRNA levels were calculated.

Validation of MXRA5 Protein
Immunocytochemistry and immunohistochemistry were used to confirm the expression of MXRA5 protein in non-fibrotic and fibrotic hCSFs and human corneas.A double immunofluorescence was performed as described earlier using the Dako mouse monoclonal anti-αSMA antibody and an Abcam goat polyclonal anti-MXRA5 antibody (ab211302) at 1:100 dilution.Invitrogen Donkey anti-Mouse IgG (H+L) Highly Cross-Adsorbed Secondary Antibody, Alexa Fluor 594 (A-21203; Thermo Fisher Scientific) and Invitrogen Donkey anti-Goat IgG (H+L) Cross-Adsorbed Secondary Antibody, Alexa Fluor 488 (A-11055; Thermo Fisher Scientific) secondary antibodies at 1:500 dilutions were used.Corneal sections were incubated in PBS for 10 minutes followed by blocking with 5% horse serum in PBST ( Jackson ImmunoResearch, West Grove, PA, USA) for 1 hour at room temperature, incubation in MXRA5 and αSMA primary antibodies overnight at 4°C, washing with PBST, and incubation with the Invitrogen Donkey anti-Mouse IgG (H+L) Highly Cross-Adsorbed Secondary Antibody, Alexa Fluor 594 (A-21203; Thermo Fisher Scientific) and Invitrogen Donkey anti-Goat IgG (H+L) Cross-Adsorbed Secondary Antibody, Alexa Fluor 488 (A-11055; Thermo Fisher Scientific) secondary antibodies for 1 hour at room temperature, and mounting with the VECTASHIELD mounting medium containing DAPI.The Leica DM4000 B fluorescence microscope equipped with a digital camera system (Spot) was used.

Statistical Analysis
Statistical analysis was performed using Prism 9.2 (Graph-Pad, Boston, MA, USA).Student's t-test was used to determine statistical significance.The data are presented as the average mean ± standard error (SEM).In statistical analysis, P ≤ 0.05 was considered significant.

Human Corneal Stromal Fibroblast and Myofibroblast Cells
The hCSFs showing the fibrotic phenotype transdifferentiated into hCMFs and acquired a myofibroblastic phenotype upon TGFβ stimulation at 72 hours.The phase-contrast microscopy and double immunofluores-FIGURE 5. Comprehensive PANTHER analysis of pathway enrichment using upregulated and downregulated genes in hCMFs compared to hCSFs.(A) The distribution of protein classes among the DEGs, distinguishing between those that were upregulated and those that were downregulated in hCMFs.(B) Biological processes used to categorize the genes highlight the variations in process participation between upregulated and downregulated genes in hCMFs.(C) Molecular function analysis of the genes at the molecular level, illustrating the functional differences in the number of genes in hCMFs.(D) Cellular component analysis showing the components in which these genes were active, highlighting the distinctions in cellular localization between upregulated and downregulated genes.Each category is presented in a format that allows for easy comparison of the functional and biological attributes associated with gene expression changes in hCMFs compared to hCSFs.cence confirmed fibroblastic features (non-fibrotic cells) (Figs.2A, 2C) and myofibroblastic features (fibrotic cells) (Figs.2B, 2D).

Qualitative and Quantitative Analyses of RNA Samples
The high purity and yield of all RNA samples extracted from hCMFs and hCSFs are shown in Table 2.The base call accuracy was 99.9% as per the Phred quality score (Q > 30) (Supplementary Fig. S1).

Differential Gene Expression Analysis
Table 3 shows the results of RNA-seq read alignment and mapping with the human reference genome (GRCh37) determined by the HISAT2.The 20,269 gene annotations were identified by StringTie and subjected to DESeq2 to obtain DEGs (hCMF vs. hCSF).The changes in gene expression between hCMFs and hCSFs indicated by the MA plot for the M values, log 2 (FC), and A values (average expressions) are shown in Supplementary Fig. S2A, and the violin plot (Supplementary Fig. S2B) unveiled differences in gene expression levels.A total of 3843 significantly expressed genes were found with P ˂ 0.05 (Supplementary Table S1).Cutoff values of log 2 (FC) ≥ 1.5 and log 2 (FC) ≤ -1.5 for upregulated and downregulated genes, respectively, led to the identification of 816 upregulated and 739 downregulated genes.The volcano plot shows that genes changed significantly (red, P ≤ 0.05), nonsignificantly (gray), and outside the cutoff values (blue, log 2 [FC]; green, P) (Fig. 3).The heatmap of the blue-yellowred color gradient (blue for highest and red for lowest) demonstrates the expression of DEGs in hCMFs and hCSFs (Fig. 4).

Pathway Enrichment Analysis
The biological relevance of 816 up-and 739 downregulated genes was scored by the PANTHER pathway enrichment classification system.Figure 5 shows the observed distributions of protein classes (Fig. 5A), variations in biological processes (Fig. 5B), differences in molecular functions (Fig. 5C), and cellular components (Fig. 5D) for the hCMFs and hCSFs.The upregulated genes in hCMFs belonged to protein classes pertaining to the cytoskeletal (42 against 18 genes), cell adhesion (29 against 12 genes), intracellular signaling proteins (41 against 28 genes), ECM (18 against seven genes), adaptor (38 against 28 genes), and translational (five against one gene).On the other hand, downregulated genes in hCMFs were in the protein classes of transcriptional regulator (51 against 76 genes), protein modifying enzyme (49 against 69 genes), and RNA metabolism (11 against 19 genes).The mapping of up-and downregulated genes was performed to understand the contributions of different biological processes in corneal microenvironment (Fig. 5B).Differences in the number of augmented proteins were found in cellular (351 against 272 genes), developmental (102 against 53 genes), metabolic (135 against 115 genes), localization (82 against 68 genes), and homeostatic (16 against 10 genes) genes, as well as suppressed proteins in the immune system (seven against 27 genes).The functional implications of DEGs in corneal wound healing were analyzed (Fig. 5C).Genes upregulated in hCMFs were linked to molecular-function activity (55 against 39), structural activity (16 against 4), transporter activity (44 against 38), and cytoskeletal motor activity (five against three).The involvement of up-and downregulated genes in TGFβinduced transdifferentiation of hCSFs to hCMFs was determined (Fig. 5D).The upregulated genes were connected to cellular anatomic entity (486 against 415 genes) and proteincontaining complex (75 against 63 genes).
The tree plot from the Reactome database exhibited hierarchical organization and connections between enriched pathways (Fig. 6).We observed changes in the interacting insulin-like growth factors, collagen assembly biosynthesis chain, activation of DNA ATR complex, cyclin A/B1/B2-associated events, and amplification from kinetochores signal genes (Fig. 6A).The dot plot revealed changes in pathways linked to the ECM, cell-cycle checkpoints, degradation of the ECM, and collagen formation (Fig. 6B).The gene clusters and combined network of comodulated genes from the enriched pathways are shown in Figure 6C.
The GSEA heatmap of the top 50 genes in hCSFs and hCMFs is presented in Figure 8A, with a red and blue color gradient (where red indicates higher and blue indicates lower expression).The most remarkably affected genes during corneal wound healing belonged to the EMT (Fig. 8B), mTORC1 signaling (Fig. 8C), angiogenesis (Fig. 8D), and TGFβ signaling pathways (Fig. 8E).

PPI Network Analysis
The PPI networks demonstrating intricate interactions and functional convergence of signaling pathways in hCSFs and hCMFs were created using Cytoscape (Supplementary Fig. S3).The maximal clique centrality (MCC) algorithm within cytoHubba, a Cytoscape plugin, examined complexity and identified the most influential nodes in PPIs.The MCC algorithm pinpointed crucial nodes within a network based on their connectivity and position within maximal cliques and revealed the nodes within the network that are important in corneal wound healing (Fig. 9A).The EMT subnetwork showed proteins such as CCN2, ACTA2, COL4A1, and TIMP3 (Fig. 9B), the mTORC1 subnetwork showed MTOR, THBS1, and SERPINE1 proteins (Fig. 9C), and the TGFβ subnetwork showed signaling proteins including TGFB1, SMAD3, SMAD7, and CTNNB1(Fig.9D).Also, the angiogenesis subnetwork exposed proteins such as PDGFA, VCAN, and SPP1, which are essential for new blood vessel formation (Fig. 9E).The integrated data from the significantly enriched pathways (EMT, mTORC1, TGFβ, and angiogenesis) are exhibited in Supplementary Table S2.Details regarding the nodes and edges of the integrated PPI network are provided in Supplementary Table S3.The top 20 influential nodes revealed by the network analysis are shown in Table 4.

MXRA Protein Expression in Human Cornea
Double immunofluorescence displaying co-localization of αSMA (red) and MXRA5 (green) in the cytoplasm of hCMFs (Figs.11F-J) but not in hCSFs (Figs.11A-E) revealed involvement of MXRA5 protein in corneal fibrosis, as well as co-localization of αSMA (red) and MXRA5 (green) protein expression in donor fibrotic human cornea (Figs.12E-H) and its absence in normal (non-fibrotic) human cornea  34 It has unique expression patterns and functional roles in maintaining the properties of human limbal stem cells. 35This action may play an essential role in improving the corneal structure and its overall function.Actin alpha-2, α-smooth muscle actin 2 (ACTA2) 6 16.63 αSMA is a hallmark of corneal wound healing and fibrosis and signifies the transdifferentiation of fibroblasts into myofibroblasts in response to TGFβ1 injury, treatment, and surgical intervention. 36Biglycan (BGN) 7 5.91 BGN is a small leucine-rich proteoglycan known to play a significant role in the structure and function of the ECM in various tissues, including the cornea.In the context of the cornea, BGN may be involved in collagen fiber assembly, contributing to the regulation of corneal transparency and biomechanical properties. 37Cellular communication network factor 2 (CCN2) 8 4.11 CCN2, or connective tissue growth factor (CTGF), promotes the proliferation and differentiation of keratocytes, the principal cells of the corneal stroma.By mediating cell adhesion, CCN2 facilitates the interaction between corneal cells and the ECM, crucial for wound healing.Additionally, its ability to enhance fibroblast growth factor-induced DNA synthesis suggests a role in corneal cell regeneration and repair following injury. 38ransforming growth factor β1 (TGFβ1) 9 4 .2 0 T G F β1 binds to its receptors (TGFBR1 and TGFBR2), initiating cellular signaling pathways that regulate various cell functions.TGFβ1 promotes keratocyte proliferation and differentiation, collagen production, and EMT processes. 39ollagen type V alpha 1 (COL5A1) 10 3.59 COL5A1 is a fibrillar-forming collagen.The interactions of type V collagen with molecules such as DNA, heparan sulfate, thrombospondin, heparin, and insulin suggest a multifaceted role in cellular functions, including cell adhesion, migration, and a possibly influence on corneal healing processes. 40hrombospondin 1 (THBS1) 11 4.93 THBS1 is an adhesive glycoprotein crucial for mediating interactions both between cells and between cells and the ECM.The role of THBS as a ligand for CD36 is particularly significant in the cornea for its antiangiogenic properties, helping maintain the avascular nature of the cornea essential for transparency and vision. 41Fibrillin 1 (FBN1) 12 3.76 FBN1 is a structural component of the ECM and forms part of the microfibrils that not only provide structural support but also play a regulatory role in the corneal tissue.These microfibrils, independent of elastin, contribute significantly to the tensile strength of the cornea and provide anchoring functions.FBN1 interacts with various growth factors, including bone morphogenetic proteins (BMPs) and latent transforming growth factor-beta-binding proteins (LTBPs), alongside cell-surface integrins and other ECM components. 42Elastin (ELN) 13 17.60 ELN is less prominent but plays a crucial role due to the unique requirements of the cornea for transparency and mechanical stability.Although elastin is a major structural protein in tissues that undergo significant expansion and recovery, its presence in the cornea is relatively minimal.However, the elastin that does exist in the cornea contributes to its ability to withstand mechanical stress and maintain its shape and structure.This is particularly important in the peripheral corneal regions and the limbus, where slight elasticity might help in accommodating intraocular pressure changes and ensuring the integrity of the ocular surface. 43ecreted protein acidic and rich in cysteine (SPARC) 14 3.28 SPARC regulates cellular behavior and maintains ECM composition through its interactions with cytokines and various components of the matrix.This multifunctional glycoprotein binds to essential molecules such as calcium, copper, diverse collagen types that are critical for corneal structure, albumin, thrombospondin (which supports cell-matrix interactions), platelet-derived growth factor (PDGF), and cell membranes.3][54] This study detected low MXRA5 mRNA and protein expression in hCSFs and normal human cornea but significantly increased levels in hCMF and fibrotic cornea (Figs.11, 12), thus revealing its role in corneal fibrosis development.Furthermore, this novel observation led to the suggestion that MXRA5 can potentially serve as an attractive target to regulate TGFβ1-induced pathological events in the cornea.The specific function of MXRA5 in corneal wound healing is still unknown.Apart from the MXRA5, increased expression of POSTN (rank 5) and COL5A1 (rank 10) genes was observed in hCMFs during RNA-seq analysis and qRT-PCR (Table 4, Fig. 10).POSTN regulates cell attachment/dispersion, fibronectin matrix, lysyl oxidase, and limbal stem cells. 34,35Such reports suggest that POSTN may play a role in maintaining corneal integrity; however, its functional role in corneal healing is still poorly known.Likewise, the functional role of COL5A1, a fibrillar-forming collagen, in corneal stroma remains elusive, despite the fact that it can interact with heparan sulfate, thrombospondin, and heparin and influence cell adhesion and migration. 40Furthermore, the current study found significant changes in WNT4, αSMA/ACTA2, AQP1, NOX4, TLL2, LMCD1, PMEPA1, LANCL2, PTX3, and many other genes (Fig. 8, Table 4) that are shown to affect TGFβ1 signaling, myofibroblast formation, reactive oxygen species production, inflammation, or ECM formation/remodeling in several non-ocular tissues [55][56][57][58][59][60] ; however, their precise role in corneal stromal healing has not yet been fully explored.
Another notable finding of this RNA-seq study is exposing a dominant role of genes associated with EMT, mTORC1 signaling, TGFβ signaling, and angiogenesis in corneal stromal fibroblast differentiation to myofibroblast (Figs.5-8).These data highlight a distinct shift in cellular mechanisms and signaling pathways during myofibroblast formation, consistent with earlier reports.Also, the PPI networks constructed from these enriched pathways show interactions and functional convergence of the pathways, serving as a testament to the complexity of corneal wound healing and highlighting the multifaceted roles of various signaling pathways in this process.
This study provides valuable new information regarding genes linked to corneal fibroblast differentiation, but it has some limitations.These include not fully addressing the impact of biological variables such as sex, age, and race; usage of human corneas from donors between 46 and 76 years old; a single dose of TGFβ1; primary cultures generated from multiple donors; and the use of 10% fetal bovine serum to obtain primary hCSF cultures.Our future studies will address these limitations by employing a larger and more diverse cohort of donors from different demographics.Also, additional research employing a broader range of TGFβ1 concentrations and serum conditions will help to better understand the spectrum of corneal stromal cellular responses in corneal wound healing and fibrosis.
In summary, this study provides a comprehensive analysis of gene count, DGE, pathways, PPI networks, and the top genes influencing corneal fibroblast differentiation to myofibroblasts.The genetic landscape governing the differences between corneal fibroblasts and myofibroblasts offers a deeper understanding of molecular dynamics of corneal wound healing and suggests potential targets for developing strategies to treat corneal fibrosis and restore vision.

FIGURE 4 .
FIGURE 4. Heatmap of gene expression differences between the hCMF and hCSF sequence profiles.This heatmap illustrates the expression levels of the genes with significant differential expression between the hCMF and hCSF sequences.The genes are displayed as rows and the hCMFs and hCSFs are shown as columns; a blue-yellow-red gradient is used to indicate the expression intensity.Blue represents higher expression, and red denotes lower expression.The heatmap also indicates key pathways associated with the DEGs.

AFIGURE 6 .
FIGURE 6. Pathway enrichment analysis using the Reactome database.(A) Tree plot of enriched pathways that illustrates the relationships between various biological pathways; the size of the nodes corresponds to the significance of the enrichment.Pathways are colored according to their enrichment scores, providing a visual representation of the key pathways involved in the DGE analysis.(B) In this dot plot, the x-axis (gene ratio) indicates the proportion of genes associated with each pathway relative to the total number of genes analyzed.The y-axis shows the names of the enriched pathways.The color gradient of the dots represents the adjusted adjusted P values, with blue indicating higher significance and red indicating lower significance.The size of the dots corresponds to the gene count involved in each pathway, with larger dots representing higher counts.(C) Network clusters showing genes involved in various ECM-related processes.The gene clusters show the genes and their interactions in five categories: ECM organization, degradation of the ECM, Assembly of collagen fibrils and other multimeric structures, non-integrin membrane-ECM interactions, and collagen biosynthesis and modifying enzymes.The integrated network combines these categories, with node size representing the number of connections (degree) and node color indicating fold change (blue for higher values, and red for lower values).

FIGURE 7 .
FIGURE 7. (A) Graphical illustration of gene set enhancement in various pathways.The most significant pathways with NES ≥ 2 and FDR q ≤ 0.05 are shown in red, and pathways with NES < 2 and/or FDR q > 0.05 are shown in green.The FDR q values of each pathway are shown in blue.(B) Heatmaps for the top four enriched pathways involving the epithelial-mesenchymal transition, mTORC1, angiogenesis, and TGFβ signaling across different samples are shown.The color gradient from blue to red indicates increasing gene expression levels.

FIGURE 8 .
FIGURE 8. Pathway enrichment analysis of hCMFs versus hCSFs.(A) Heatmap of the top 50 genes in each cell line, with the red-to-blue color scale indicating higher to lower gene expression.(B-E) Significantly enriched pathways: EMT pathway (NES = 5.569), mTORC1 signaling pathway (NES = 2.949), angiogenesis (NES = 2.176), and TGFβ signaling (NES = 2.008).These pathways highlight a shift toward increased cellular differentiation and fibrotic responses in hCMFs, suggesting potential impacts on corneal fibrosis and transparency.

FIGURE 9 .
FIGURE 9. PPI networks.(A) Integrated protein-protein interaction network of significantly enriched pathways in hCMFs.(B-E) Subnetworks of the four pathways related to the genes exhibiting the most significant enrichment: EMT (B), mTORC1 signaling (C), TGFβ signaling (D), and angiogenesis (E), derived from GSEA using gene counts.Influential nodes, as determined by the MCC algorithm using CytoHubba in Cytoscape, are emphasized to indicate their central role in the network (red to yellow color scale, with red indicating the most influential nodes).

FIGURE 11 .
FIGURE 11.Double immunofluorescence staining results for αSMA and MXRA5 in non-fibrotic and fibrotic corneal stromal cells.(A-E) Non-fibrotic cells.(F-J) Fibrotic cells.Panel A shows a differential interference contrast (DIC) image depicting the general morphology of non-fibrotic cells, and panel F shows fibrotic cells.Panels B and G display the immunofluorescence images for αSMA (red), which is absent in non-fibrotic cells but clearly expressed in fibrotic cells.Panels C and H show the immunofluorescence images for MXRA5 (green), indicating its expression in both cell types but with higher expression in fibrotic cells.Panels D and I display DAPI staining (blue) marking the nuclei of the cells.Panels E and J are merged images combining αSMA, MXRA5, and DAPI.The comparison highlights the elevated expression of MXRA5 in fibrotic cells, as indicated by the more intense green fluorescence in panels H and J. (K) The quantification revealed significant (P < 0.001) increase in expression of MXRA5 and αSMA in corneal stromal fibrotic cells compared to non-fibrotic cells.Scale bar: 40 μm.

TABLE 2 .
RNA Sample Preparation Details for hCSF and hCMF Samples Using the RNeasy Mini Kit

TABLE 3 .
Mapping Summary of the Eight Sequence Samples

TABLE 4 .
Influential Nodes With a Possible Role in Corneal Wounds and Their Ranking in the Network 44

TABLE 4 .
51ntinued key role in regulating the formation of collagen fibrils.As a member of the small leucine-rich proteoglycan (SLRP) family, specifically within the SLRP class II subfamily, it significantly contributes to collagen fibrillogenesis, impacting the mechanical properties and overall health of the cornea.51

TABLE 5 .
Tope Ten Upregulated and Downregulated Differentially Expressed Genes the best of our knowledge, to the best of our knowledge, identified expression of the MXRA5 gene in human corneal stroma and its role in corneal fibrosis and the involvement of POSTN and COL5A1 genes in human corneal stromal cell modulation.The non-ocular literature reveals the expression and functional role of MXRA5 in TGFβ1-driven fibrotic processes